home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / basic / _matrix.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  7.4 KB  |  416 lines

  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  _matrix.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. *******************************************************************************/
  13.  
  14.  
  15.  
  16. #include <LEDA/matrix.h>
  17.  
  18. #include <math.h>
  19.  
  20. #define EPSILON 1e-12
  21.  
  22. //------------------------------------------------------------------------------
  23. //  matrix member functions
  24. //------------------------------------------------------------------------------
  25.  
  26.  
  27. void matrix::flip_rows(int i,int j)
  28. { vector* p = v[i];
  29.   v[i] = v[j];
  30.   v[j] = p;
  31.  }
  32.  
  33.  
  34. matrix::matrix(int dim1, int dim2)  
  35. {
  36.   if (dim1<0 || dim2<0) 
  37.   error_handler(1,"matrix: negative dimension."); 
  38.  
  39.   d1=dim1; 
  40.   d2=dim2; 
  41.  
  42.   if (d1 > 0) 
  43.   { v = new vector*[d1];
  44.     for (int i=0;i<d1;i++) v[i] = new vector(d2); 
  45.    }
  46.   else v = nil;
  47. }
  48.  
  49.  
  50. matrix::matrix(const matrix& p)  
  51.   d1 = p.d1;
  52.   d2 = p.d2;
  53.     
  54.   if (d1 > 0) 
  55.   { v = new vector*[d1];
  56.     for (int i=0;i<d1;i++) v[i] = new vector(*p.v[i]); 
  57.    }
  58.   else v = nil;
  59. }
  60.  
  61. matrix::matrix(int dim1, int dim2, double** p)  
  62. { d1=dim1; 
  63.   d2=dim2; 
  64.   v = new vector*[d1];
  65.   for(int i=0;i<d1;i++) 
  66.   { v[i] = new vector(d2); 
  67.     for(int j=0;j<d2;j++) elem(i,j) = p[i][j];
  68.    }
  69.  
  70.  }
  71.  
  72. matrix::~matrix()  
  73. { if (v) 
  74.   { while(d1--) delete v[d1]; 
  75.     delete v;
  76.    }
  77. }
  78.  
  79.  
  80. void matrix::check_dimensions(const matrix& mat) const
  81. { if (d1 != mat.d1 || d2 != mat.d2)
  82.    error_handler(1,"incompatible matrix types.");
  83.  }
  84.  
  85. matrix::matrix(const vector& vec)
  86. { d1 = vec.d;
  87.   d2 = 1;
  88.   v = new vector*[d1];
  89.   for(int i=0; i<d1; i++)
  90.   { v[i] = new vector(1);
  91.     elem(i,0) = vec[i];
  92.    }
  93.     
  94. }
  95.  
  96. matrix& matrix::operator=(const matrix& mat)
  97. { register int i,j;
  98.  
  99.   if (d1 != mat.d1 || d2 != mat.d2)
  100.   { for(i=0;i<d1;i++) delete v[i];
  101.     delete v;
  102.     d1 = mat.d1;
  103.     d2 = mat.d2;
  104.     v = new vector*[d1];
  105.     for(i=0;i<d1;i++) v[i] = new vector(d2);
  106.    }
  107.  
  108.   for(i=0;i<d1;i++)
  109.     for(j=0;j<d2;j++) elem(i,j) = mat.elem(i,j);
  110.  
  111.   return *this;
  112. }
  113.  
  114. int matrix::operator==(const matrix& x) const
  115. { register int i,j;
  116.   if (d1 != x.d1 || d2 != x.d2) return false;
  117.  
  118.   for(i=0;i<d1;i++)
  119.     for(j=0;j<d2;j++)
  120.       if (elem(i,j) != x.elem(i,j)) return false;
  121.  
  122.   return true;
  123.  }
  124.  
  125.  
  126. vector& matrix::row(int i) const
  127. { if ( i<0 || i>=d1 )  error_handler(1,"matrix: row index out of range");
  128.    return *v[i];
  129. }
  130.  
  131.  
  132. double& matrix::operator()(int i, int j)
  133. { if ( i<0 || i>=d1 )  error_handler(1,"matrix: row index out of range");
  134.   if ( j<0 || j>=d2 )  error_handler(1,"matrix: col index out of range");
  135.   return elem(i,j);
  136. }
  137.  
  138. double matrix::operator()(int i, int j) const
  139. { if ( i<0 || i>=d1 )  error_handler(1,"matrix: row index out of range");
  140.   if ( j<0 || j>=d2 )  error_handler(1,"matrix: col index out of range");
  141.   return elem(i,j);
  142. }
  143.  
  144. vector matrix::col(int i)  const
  145. { if ( i<0 || i>=d2 )  error_handler(1,"matrix: col index out of range");
  146.   vector result(d1);
  147.   int j = d1;
  148.   while (j--) result.v[j] = elem(j,i);
  149.   return result;
  150. }
  151.  
  152. matrix::operator vector() const
  153. { if (d2!=1) 
  154.    error_handler(1,"error: cannot make vector from matrix\n");
  155.   return col(0);
  156. }
  157.  
  158. matrix matrix::operator+(const matrix& mat)
  159. { register int i,j;
  160.   check_dimensions(mat);
  161.   matrix result(d1,d2);
  162.   for(i=0;i<d1;i++)
  163.    for(j=0;j<d2;j++)
  164.       result.elem(i,j) = elem(i,j) + mat.elem(i,j);
  165.   return result;
  166. }
  167.  
  168. matrix& matrix::operator+=(const matrix& mat) 
  169. { register int i,j;
  170.   check_dimensions(mat);
  171.   for(i=0;i<d1;i++)
  172.    for(j=0;j<d2;j++)
  173.       elem(i,j) += mat.elem(i,j);
  174.   return *this;
  175. }
  176.  
  177. matrix& matrix::operator-=(const matrix& mat) 
  178. { register int i,j;
  179.   check_dimensions(mat);
  180.   for(i=0;i<d1;i++)
  181.    for(j=0;j<d2;j++)
  182.       elem(i,j) -= mat.elem(i,j);
  183.   return *this;
  184. }
  185.  
  186.  
  187. matrix matrix::operator-(const matrix& mat)
  188. { register int i,j;
  189.   check_dimensions(mat);
  190.   matrix result(d1,d2);
  191.   for(i=0;i<d1;i++)
  192.    for(j=0;j<d2;j++)
  193.       result.elem(i,j) = elem(i,j) - mat.elem(i,j);
  194.   return result;
  195. }
  196.  
  197.  
  198. matrix matrix::operator-()  // unary
  199. { register int i,j;
  200.   matrix result(d1,d2);
  201.   for(i=0;i<d1;i++)
  202.    for(j=0;j<d2;j++)
  203.       result.elem(i,j) = -elem(i,j);
  204.   return result;
  205. }
  206.  
  207.  
  208. matrix matrix::operator*(double f)
  209. { register int i,j;
  210.   matrix result(d1,d2);
  211.   for(i=0;i<d1;i++)
  212.    for(j=0;j<d2;j++)
  213.       result.elem(i,j) = elem(i,j) *f;
  214.   return result;
  215. }
  216.  
  217. matrix matrix::operator*(const matrix& mat)
  218. { if (d2!=mat.d1)
  219.      error_handler(1,"matrix multiplication: incompatible matrix types\n");
  220.   
  221.   matrix result(d1, mat.d2);
  222.   register int i,j;
  223.  
  224.   for (i=0;i<mat.d2;i++)
  225.   for (j=0;j<d1;j++) result.elem(j,i) = *v[j] * mat.col(i);
  226.  
  227.  return result;
  228.  
  229. }
  230.  
  231. double matrix::det() const
  232. {
  233.  if (d1!=d2)  
  234.    error_handler(1,"matrix::det: matrix not quadratic.\n");
  235.  
  236.  int n = d1;
  237.  
  238.  matrix M(n,1);
  239.  
  240.  int flips;
  241.  
  242.  double** A = triang(M,flips);
  243.  
  244.  if (A == NULL)  return 0;
  245.  
  246.  double Det = 1;
  247.  
  248.  for(int i=0;i<n;i++) Det *= A[i][i];
  249.  
  250.  for(i=0;i<n;i++) delete A[i];
  251.  delete A;
  252.  
  253.  return (flips % 2) ? -Det : Det;
  254.  
  255. }
  256.  
  257.  
  258. double** matrix::triang(const matrix& M, int& flips)  const
  259. {
  260.  register double **p, **q;
  261.  register double *l, *r, *s;
  262.  
  263.  register double pivot_el,tmp;
  264.  
  265.  register int i,j, col, row;
  266.  
  267.  int n = d1;
  268.  int d = M.d2;
  269.  int m = n+d;
  270.  
  271.  double** A = new double*[n];
  272.  
  273.  p = A;
  274.  
  275.  for(i=0;i<n;i++) 
  276.  { *p = new double[m];
  277.    l = *p++;
  278.    for(j=0;j<n;j++) *l++ = elem(i,j);
  279.    for(j=0;j<d;j++) *l++ = M.elem(i,j);
  280.   }
  281.  
  282.  flips = 0;
  283.  
  284.  for (col=0, row=0; row<n; row++, col++)
  285.  { 
  286.    // search for row j with maximal absolute entry in current col
  287.    j = row;
  288.    for (i=row+1; i<n; i++)
  289.      if (fabs(A[j][col]) < fabs(A[i][col])) j = i;
  290.  
  291.    if ( n > j && j > row) 
  292.    { double* p = A[j];
  293.      A[j] = A[row];
  294.      A[row] = p;
  295.      flips++;
  296.     }
  297.  
  298.    tmp = A[row][col];
  299.    q  = &A[row];
  300.  
  301.    if (fabs(tmp) < EPSILON) // matrix has not full rank
  302.    { p = A;
  303.      for(i=0;i<n;i++) delete *p;
  304.      delete A;
  305.      return NULL;
  306.     }
  307.  
  308.    for (p = &A[n-1]; p != q; p--)
  309.    { 
  310.      l = *p+col;
  311.      s = *p+m;    
  312.      r = *q+col;
  313.  
  314.      if (*l != 0.0)
  315.      { pivot_el = *l/tmp;
  316.         while(l < s) *l++ -= *r++ * pivot_el;
  317.       }
  318.  
  319.     }
  320.  
  321.   }
  322.  
  323.  return A;
  324. }
  325.  
  326. matrix matrix::inv() const
  327. {
  328.  if (d1!=d2)  
  329.      error_handler(1,"matrix::inv: matrix not quadratic.\n");
  330.  int n = d1;
  331.  matrix I(n,n);
  332.  for(int i=0; i<n; i++) I(i,i) = 1;
  333.  return solve(I);
  334. }
  335.  
  336.  
  337.  
  338. matrix matrix::solve(const matrix& M) const
  339. {
  340.  
  341. if (d1 != d2 || d1 != M.d1)
  342.      error_handler(1, "Solve: wrong dimensions\n");
  343.  
  344.  register double **p, ** q;
  345.  register double *l, *r, *s;
  346.  
  347.  int      n = d1;
  348.  int      d = M.d2;
  349.  int      m = n+d;
  350.  int      row, col,i;
  351.  
  352.  
  353.  double** A = triang(M,i);
  354.  
  355.  if (A == NULL) 
  356.    error_handler(1,"matrix::solve: matrix has not full rank.");
  357.  
  358.  for (col = n-1, p = &A[n-1]; col>=0; p--, col--)
  359.  { 
  360.    s = *p+m;
  361.  
  362.    double tmp = (*p)[col];
  363.  
  364.    for(l=*p+n; l < s; l++) *l /=tmp;
  365.  
  366.    for(q = A; q != p; q++ )
  367.    { tmp = (*q)[col];
  368.      l = *q+n;
  369.      r = *p+n;
  370.      while(r < s)  *l++ -= *r++ * tmp;
  371.     }
  372.                  
  373.   } 
  374.  
  375.   matrix result(n,d);
  376.  
  377.   for(row=0; row<n; row++)
  378.   { l = A[row]+n;
  379.     for(col=0; col<d; col++) result.elem(row,col) = *l++;
  380.     delete A[row];
  381.    }
  382.  
  383.   delete A;
  384.  
  385.   return result;
  386. }
  387.  
  388.  
  389.  
  390.  
  391. matrix matrix::trans() const
  392. { matrix result(d2,d1);
  393.   for(int i = 0; i < d2; i++)
  394.     for(int j = 0; j < d1; j++)
  395.       result.elem(i,j) = elem(j,i);
  396.   return result;
  397. }
  398.      
  399.  
  400. ostream& operator<<(ostream& s, const matrix& M)
  401. { int i;
  402.   s <<"\n";
  403.   for (i=0;i<M.d1;i++) s << M[i] << "\n"; 
  404.   return s;
  405. }
  406.  
  407. istream& operator>>(istream& s, matrix& M)
  408. { int i=0;
  409.   while (i<M.d1 && s >> M[i++]);
  410.   return s;
  411. }
  412.  
  413.  
  414.